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ABSTRACT 

A  discontinuous  Galerkin  formulation  that  avoids  the  use  of  discrete  quadrature  formulas  is 
described  and  applied  to  linear  and  nonlinear  test  problems  in  one  and  two  space  dimensions.  This 
approach  requires  less  computational  time  and  storage  than  conventional  implementations  but 
preserves  the  compactness  and  robustness  inherent  in  the  discontinuous  Galerkin  method.  Test 
problems  include  the  linear  and  nonlinear  one-dimensional  scalar  advection  of  both  smooth  and 
discontinuous  initial  value  problems,  the  two-dimensional  scalar  advection  of  smooth  initial  value 
problems  that  are  discretized  by  using  unstructured  grids  with  varying  degrees  of  smoothness 
and  regularity,  and  two-dimensional  linear  Euler  solutions  on  unstructured  grids. 
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Introduction 


Computational  methods  for  aeroacoustics  must  possess  accuracy  properties  that  exceed 
those  of  conventional  second-order  computational  fluid  dynamics  (CFD)  methods.  At  the  same 
time,  many  problems  of  interest  involve  complex  geometries  that  are  not  easily  treated  by 
common  high-order  methods  that  normally  require  a  smooth,  structured  grid.  In  addition  to  the 
geometrically  complex  problem,  we  are  particularly  interested  in  strongly  nonlinear  flows  that 
contain  shock  waves  as  a  major  source  of  sound  generation,  such  as  in  the  case  of  jet  noise. 

In  an  effort  to  satisfy  these  requirements,  the  relatively  untried  discontinuous  Galerkin  (DG) 
method  is  being  tested  for  hyperbolic  problems.  Some  advantages  of  this  approach  include 
the  ease  with  which  the  method  can  be  applied  to  both  structured  and  unstructured  grids  and  its 
suitability  for  parallel  computer  architectures.  This  approach  also  has  several  useful  mathematical 
properties. 

Johnson  and  PitkSrata^  proved  stability  and  error  estimates  for  linear  scalar  advection.  In 
a  series  of  papers,  Cockbum  et  al.^"^  discussed  the  DG  method  with  the  use  of  approximate 
Riemann  solvers,  limiters,  and  total-variation  diminishing  (TVD)  Runge-Kutta  time  discretiza¬ 
tions  for  nonlinear  hyperbolic  problems.  In  reference  2,  the  general  formulation  in  one  space 
dimension  is  provided;  This  work  includes  an  analysis  of  the  accuracy  and  stability  (in  terms 
of  total  variation)  and  a  detailed  description  of  the  implementation.  Numerical  examples  using 
the  scalar  equations  are  also  provided.  In  reference  3,  the  method  is  applied  to  one-dimensional 
systems.  Stability  of  the  initial  boundary  value  problem  for  a  linear  system  is  proved,  and  nu¬ 
merical  experiments  are  presented  for  the  one-dimensional  Euler  equations.  In  reference  4,  the 
method  is  generalized  to  multispace  dimensions.  A  main  result  in  reference  4  is  the  design  of 
a  limiter  that  applies  to  general  triangulations,  maintains  a  high  order  of  accuracy  in  smooth 
regions,  and  guarantees  maximum  norm  stability.  Jiang  and  Shu^  have  also  proved  that  the  DG 
method  satisfies  a  local  cell  entropy  inequality  for  the  square  entropy,  for  arbitrary  triangulations 
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in  any  space  dimension,  and  for  any  order  of  accuracy.  This  proof  trivially  implies  stability 
of  the  method  for  nonlinear  shocked  problems  in  the  scalar  case. 

Although  the  DG  method  has  not  been  widely  used  in  the  CFD  arena,  several  instances 
exist  in  which  this  method  has  been  applied  to  the  Euler  or  Navier-Stokes  equations.^’"^’^^  Halt 
and  Agarwal®  applied  the  method  of  moments  (similar  to  the  DG  method)  to  the  steady,  two- 
dimensional  Euler  equations  for  subsonic  flows.  Bassi  and  Rebay^  applied  the  DG  method 
to  two-dimensional  Euler  equations  for  transonic  flows  and  demonstrated  the  importance  of  the 
proper  treatment  of  curved  boundaries.  In  reference  8,  Bassi  and  Rebay  extended  their  method  to 
the  Navier-Stokes  equations  by  introducing  the  gradient  of  the  solution  as  an  auxiliary  variable. 
Most  recently,  Lowrie,  Roe,  and  van  Leer^  presented  a  fully  discrete  DG  method  for  the  unsteady 
Euler  equations.  Biswas,  Devine,  and  Flaherty'®  applied  the  DG  method  in  an  h-p  version  in 
an  adaptive  grid  environment  and  considered  the  issues  in  using  limiters  for  moments,  as  well 
as  for  parallel  implementations. 

In  the  works  described  above,  the  integrals  that  appear  in  the  formulation  are  evaluated  with 
quadrature  formulas.  In  this  particular  work,  the  DG  method  is  implemented  in  a  form  that 
avoids  the  use  of  quadrature  formulas  with  the  goal  of  improving  the  efficiency  of  the  method. 
In  the  first  section,  the  DG  method  is  described  in  general.  Then,  additional  motivation  for  using 
the  quadrature-free  approach  is  given,  along  with  specific  details  of  the  formulation.  Storage 
requirements  and  operation  count  are  also  discussed,  and  the  results  of  a  stability  analysis  are 
given.  It  is  shown  that  the  primary  arguments  against  the  method,  high  storage  requirements  and 
a  high  operation  count,  are  unjustified.  The  last  section  presents  numerical  results  for  one-  and 
two-dimensional  test  problems.  A  linear  scalar  equation  is  used  to  verify  the  general  properties 
of  the  DG  method  for  solution  expansions  up  to  degree  12.  The  nonlinear  Burgers’  equation  is 
used  to  demonstrate  the  shock-capturing  capabilities  of  the  method  in  one  space  dimension.  The 
method  is  applied  to  both  scalar  advection  and  the  linear  Euler  equations  in  two  space  dimensions 
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to  demonstrate  the  unstructured  grid  capability.  This  work  focuses  on  the  new  formulation  of 
the  DG  method  and  defers  all  discussion  of  boundary  conditions  to  future  articles.  Hence,  all 
test  problems  are  treated  as  spatially  periodic. 

General  Discontinuous  Galerkin  Method 

Consider  an  arbitrary  domain  fl  in  which  the  solution  is  governed  by  a  conservation  equation 
of  the  form 

Ut  +  V  ■  F  =  0  (1) 

The  DG  method  can  be  obtained  by  partitioning  the  domain  onto  smaller,  nonoverlapping 
elements  0.i  that  cover  the  domain  and  then  applying  a  traditional  Galerkin**  method  to  each 
element.  The  Galerkin  approach  within  an  element  is  defined  by  selecting  a  finite-dimensional 
basis  set,  approximating  the  solution  as  an  expansion  in  that  basis,  and  then  projecting  the 
governing  equation  onto  each  member  of  that  basis  set: 

B  =  {bk,  0  <  k  <  N{p,d)}  (2) 

N 

^  ^  ^  '^i,j  bj  (3) 

J  bk(Vt  +  V  ■  pydn  =  0  {0  <  k  <  N,  0  <  i  <  I)  (4) 

Qi 

where  I  is  the  total  number  of  elements.  The  basis  set  must  be  constructed  from  the  lower 
order  terms  of  a  complete  and  linearly  independent  set.  In  the  fully  discrete  approach,^  the  basis 
set  contains  both  temporal  and  spatial  functions.  In  the  semidiscrete  approach,  which  is  used 
here,  the  basis  set  contains  only  spatial  functions,  and  the  solution  expansion  coefficients  vij  are 
functions  of  time.  The  number  of  terms  in  the  expansion,  N{p,d)  -|- 1,  depends  on  the  degree 
of  the  expansion  p  and  the  number  of  time-space  dimensions  d  that  are  represented.  When  the 
basis  functions  are  polynomials  of  degree  p,  Johnson  et  al.*  proved  the  order  of  accuracy  of  the 
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method  to  be  at  least  p+1/2.  In  practice,  however,  the  order  of  accuracy  of  the  method  observed 
in  most  cases  is  p  + 1.  Except  where  the  full  form  is  needed  for  clarity,  N{p,  d)  will  be  denoted 
simply  by  N  to  streamline  the  notation. 

The  divergence  term  is  recast  by  applying  integration  by  parts  as 

j  hY,{vi^,\bjd(l  -  J  Vbk  ■  F^dn  +  J  bkF,^  ■  ds  =  0 

fii  dQi 

{0  <  k  <  N,  0  <  i  <  I)  (5) 


where  ds  is  an  outward-pointing  surface-element  normal  and  the  flux  vector  F^  is  some 
approximation  to  a  Riemann  flux  that  depends  on  the  solutions  in  both  element  i  and  the 
neighboring  element.  The  approximate  Riemann  flux  provides  the  crucial  coupling  between 
elements  as  well  as  the  correct  upwind  bias  required  to  ensure  stability. 

Although  equation  (5)  could  be  evaluated  in  physical  space,  the  equation  can  be  conveniently 
(and  sometimes  advantageously)  represented  in  terms  of  coordinates  that  are  local  to  the  element. 
After  such  a  mapping,  equation  (5)  becomes 

J  bkY^{v^,J\bjJidn  -  js7bk-  iJ^FJidQ  +  J  bkJ-^Fi^  ■  J^ds  =  0 

Qi  Qi  dUi 

(0  <  A:  <  iV,  0  <  z  <  /) 

(6) 


where 


d{x,y,z) 


Ji 


|Ji 


*  "  ’ 

and  (^,  p,  C)  are  the  local  coordinates  within  the  element. 

Depending  on  the  choice  of  basis  functions  and  the  form  of  Ji,  the  first  integral  can  usually 
be  simplified.  If,  for  example,  the  basis  set  is  orthonormal  with  respect  to  Ji,  then  the  first 
integral  term  in  equation  (6)  reduces  to  an  identity  mass  matrix.  In  most  implementations,  either 
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the  transformation  is  an  isoparametric  form  (thus,  Ji  is  a  polynomial)  or  the  element  is  allowed 
to  be  an  arbitrary  polygon.  In  either  case,  the  temporal  term  can  be  expressed  as  M,  V* ,  where 
Mi  =  \rnk,j\  is  the  mass  matrix,  V  =  [u/],  and 

=  /  bkbjJidn  (7) 

In  most  cases,  the  mass  matrix  can  be  computed  and  inverted  in  advance  of  the  main  calculation; 
however,  it  must  be  stored  for  every  element.  In  all  previous  implementations,  the  spatial 
integrals  are  evaluated  with  quadrature  formulas  that  are  appropriate  to  the  element  shape  and 
the  required  degree  of  accuracy.  The  resulting  +  1  equations  are  then  used  to  temporally 
evolve  the  coefficients  of  the  solution  expansion  vij. 


In  the  present  formulation,  we  place  a  fundamental  restriction  on  the  types  of  elements  that 
are  permitted.  In  particular,  we  require  that  all  elements  (except  perhaps  elements  near  a  curved 
boundary)  have  a  linear  mapping  to  a  simple  similarity  element,  such  as  an  equilateral  triangle  or 
a  square  in  two  dimensions  (or  other  simple  elements  in  three  dimensions).  With  this  restriction, 
the  temporal  term  can  be  rewritten  as  JjM  ,  where  the  mass  matrix  M  is  now  defined  by 


mk,j  = 


(8) 


As  a  consequence,  M  is  the  same  for  all  elements  of  a  given  similarity  shape  (e.g.,  one  M  is 
applicable  to  all  triangles,  another  is  applicable  to  all  squares,  and  so  on)  except  those  near  a 
curved  boundary.  The  considerable  reduction  in  the  amount  of  storage  needed  far  outweighs  the 
inconvenience  that  may  arise  from  the  restriction  in  element  shape. 

Throughout  the  remainder  of  this  discussion,  the  subscript  that  identifies  the  element  will 
be  omitted  unless  its  use  is  necessary  for  clarity.  In  addition,  although  the  above  discussion 
applies  to  a  general  basis  set,  the  remaining  discussion  will  assume  that  the  basis  functions  are 
polynomials.  In  particular,  the  results  presented  later  use  a  basis  that  is  constructed  from  simple 
monomials  (e.g.,  B  =  {l,  x,  y,  xy,  y^^...}). 
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Quadrature-Free  Approach 


Quadrature  formulas,  especially  Gaussian  quadrature  formulas,  are  usually  the  most  accurate 
and  efficient  means  of  evaluating  integrals.  However,  this  feature  is  based  on  the  assumption 
that  the  data  are  readily  available  at  the  quadrature  points  (i.e.,  the  data  are  stored  there),  which 
is  typical  for  most  finite-element  and  spectral-element  methods.  In  these  methods,  the  unknown 
variables  are  the  values  of  the  solution  at  the  quadrature  points,  and  the  basis  or  “shape”  functions 
are  often  the  Lagrangian  polynomials  associated  with  the  quadrature  points.  In  the  DG  method, 
however,  the  formulation  requires  the  evaluation  of  both  volume  and  surface  integrals,  and  no 
single  set  of  iV  -I- 1  quadrature  points  exists  that  can  be  used  to  evaluate  all  integrals  to  the 
required  accuracy.  Thus,  the  usual  practice  in  the  implementation  of  the  DG  method  is  to 
choose  a  basis  and  store  the  expansion  coefficients.  As  a  consequence,  the  evaluation  of  the 
volume  integral,  for  example,  requires  N  +  1  operations  at  each  quadrature  point  simply  to 
obtain  the  data  needed  to  evaluate  the  quadrature  formula.  Because  all  integrals  must  be  exact 
for  polynomials  of  degree  2p,  the  operation  count  for  the  complete  evaluation  of  the  volume 
integral  is  on  the  order  of  2  to  3  times  A''  -1- 1  operations  for  each  equation,  with  the  assumption 
that  optimal  quadrature  formulas  exist. 

In  contrast,  if  the  basis  functions  are  judiciously  chosen  such  that  integrals  of  products  of 
the  basis  functions  can  be  evaluated  exactly,  then  the  complete  volume  integral  can  be  evaluated 
exactly  in  iV  -f  1  or  fewer  operations.  A  more  precise  comparison  of  the  operation  count  is 
given  in  a  later  section.  In  the  quadrature-free  implementation,  we  derive  a  set  of  matrices 
that  accomplish  this  task.  Further,  as  a  consequence  of  the  restriction  previously  placed  on 
the  element  type,  these  matrices  are  the  same  for  all  elements  of  a  given  type  (i.e.,  one  set  of 
matrices  applies  to  all  triangles,  another  set  applies  to  all  squares,  and  so  on).  The  result  is  a 
method  with  low  storage  requirements  and  low  computational  cost  that  retains  the  ability  to  treat 
unstructured  grids  in  an  accurate  and  robust  manner. 
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Flux  Expansion 


To  accomplish  this  task,  however,  an  efficient  method  for  expanding  the  flux  vector  F  in 
terms  of  the  basis  set  is  needed: 

M 

F  =  M  >  N{p,d)  (9) 

j=0 

where  the  number  of  terms  in  the  expansion  M  depends  on  the  form  of  the  nonlinearity  in  F. 
To  obtain  the  design  order  of  accuracy,  M  must  be  at  least  as  large  as  N{p+l,d),  such  that 
the  degree  of  v6jt  •  F  is  at  least  as  large  as  the  degree  of  bkV. 

The  expansion  is  trivially  accomplished  when  F  is  a  linear  function  of  U.  In  this  case,  the 
coefficients  fj  are  identically  zero  for  j  >  N{p,d).  Similarly,  for  common  test  problems  such 
as  the  nonlinear  Burgers’  equation,  F  =  ^U[x)  can  be  obtained  by  multiplication  so  long  as 
triple  (and  higher)  products  of  the  basis  function  are  also  easily  integrated.  Such  is  the  case 
when  the  basis  functions  are  polynomials  in  the  local  coordinates  of  a  simple  element. 

The  more  complex  flux  functions,  such  as  those  of  the  Euler  equations,  can  be  treated 
in  several  ways.  One  approach  is  that  of  Lowrie  et  al.,^  in  which  the  solution  expansion  is 
applied  to  the  parameter  vector  y/pu^  y/pv,  y/pH^]  rather  than  the  conserved  variables 
[p,  pu,  pv,  pe].  Although  this  approach  allows  all  nonlinear  fluxes  to  be  evaluated  exactly, 
the  coefficient  matrices  become  dependent  on  the  local  solution.  Hence,  the  coefficient  matrices 
derived  here  would  be  different  for  every  element  and  would  need  to  be  recomputed  at  every  time 
step;  as  a  result,  the  storage  requirements  and  computational  effort  would  be  unacceptably  high. 

A  general  alternative  is  to  expand  the  flux  in  a  Taylor  series.  Because  the  Euler  equations 
consist  of  terms  such  as  (pu)  j p  and  {pu){pv)l p,  a  practical  procedure  is  to  approximate  p~ 
in  the  basis  functions  by  using  a  Taylor  series  and  then  use  multiplication  to  complete  the  flux 
expansion. 
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Another  approach  is  to  define  the  flux  in  terms  of  the  projection  operator  as  follows.  Let 

M 

^  p  =  rjbj  (10) 

j-0 


then 


P  V 


=  1 


J  bkppdCl  =  J  bkdO,  (k  =  0,1,  2,  — ,  M) 

n  Q 


(11) 


This  set  of  equations  results  in  a  system  that  is  linear  in  rj  and  can  easily  be  solved  exactly. 
The  flux  terms,  such  as  {pu){pv)lp,  can  then  be  computed  by  multiplication,  just  as  in  the 
Taylor  series  approach.  Alternatively,  the  projection  method  can  be  used  to  determine  the  flux 
expansion  directly.  Let 


Then, 


M 

{pu){pv)lp  «  F  =  ^bjfj 
i=o 


(12) 


JbkpFdn  =  J bk{pu){pv)dn  {k  =  0,1,2,...,  M)  (13) 

n  SI 


defines  a  set  of  equations  that  is  linear  for  fj. 

The  projection  approach  is  appealing  because  the  error  of  the  approximate  flux  is  minimized 
over  the  whole  of  the  element,  whereas  in  the  Taylor  series  approach  the  flux  will  be  more 
accurate  near  the  center  of  the  element  than  near  the  edges  of  the  element. 

Volume  Integral 


After  the  flux  is  represented  as  in  equation  (9),  the  volume  integral  of  equation  (6)  can 
easily  be  rewritten  as  a  matrix  times  a  vector  as 

JMVf  -  G  -  F  +  j  •  Jds  =  0 


da 


(14) 


where  V  =  [vj.],  G 


,  B  = 

9k,i  =  /  =  jj-'f, 

o 

II 

o 

II 

(15) 
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and  M  >  N  depends  on  the  degree  to  which  a  nonlinear  flux  must  be  expanded.  The  vector 
matrix  G  is  a  constant  for  all  elements  of  a  given  family  and  can  be  precomputed  and  stored. 
The  evaluation  of  the  volume  integral  becomes  an  order  M  +  l  operation  for  each  of  the  iV  +  1 
equations  in  the  element. 

Boundary  Integral 

The  boundary  integral  is  partitioned  into  segments  associated  with  the  sides  of  the  element. 
The  integral  on  each  boundary  segment  can  be  rewritten  in  terms  of  a  vector  times  a  precomputed 
matrix;  however,  the  procedure  is  complicated  by  the  fact  that  is  a  function  of  the  solution 
in  the  two  elements  on  either  side  of  a  boundary  segment  and  each  of  these  elements  has  a 
distinct  local  coordinate  system.  The  remedy  is  to  translate  each  basis  function  from  its  local 
coordinate  system  to  a  coordinate  system  that  is  common  to  both  elements.  The  use  of  a 
boundary-segment-based  coordinate  system  actually  results  in  a  reduction  in  the  total  work  and 
storage  required.  The  complete  process  can  be  summarized  in  three  steps:  translate  the  solution 
to  an  edge  coordinate  system,  compute  the  approximate  Riemann  flux  in  the  edge  coordinates, 
and  project  the  flux  onto  the  space  defined  by  the  element  basis  set. 

For  clarity,  the  procedure  is  illustrated  for  a  triangular  element;  however,  the  same  procedure 
can  be  applied  to  all  element  types,  and  the  procedure  facilitates  the  use  of  mixed  element  types 
(i.e.,  squares  and  triangles  together).  Figure  1  illustrates  a  general  triangle  that  has  been  mapped 
into  a  equilateral  triangle.  The  equilateral  triangle  has  a  local  coordinate  system  rj)  with  its 
origin  at  the  centroid  of  the  element;  each  edge  also  has  a  local  coordinate  ^  with  its  origin  in 
the  center  of  the  edge.  Note  that  the  dimensionality  of  the  edge  coordinate  system  is  lower 
than  that  of  the  element  by  1. 

For  each  edge  of  the  triangle,  a  constant  matrix  Tj  exists  that  relates  each  member  of  the 
local  basis  to  an  expansion  in  terms  of  the  edge  coordinate.  The  subscript  j  identifies  the  edge 
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X  —  Xo 

.y-Vo 


=  J 


Figure  1.  Transform  from  general  triangle  to  equilateral  similarity  triangle  that  shows 

coordinate  systems  associated  with  interior  and  edges  of  similarity  triangle, 
to  which  the  matrix  applies.  For  example,  on  edge  0, 


B  =  [bi] 


■1  - 

e 

V 

e 

2 


1 

0 

w 

0 

i 

12 


0 

1 

0 

0 

-1 

2V3 

0 


0 

0 

0 

1 

0 

0 


=  To  bj  =  To  B 


(16) 


The  matrices  for  the  other  two  edges  are  considerably  more  dense  but  are  easily  derived  in  exact 
form  with  the  aid  of  a  symbolic  algebra  package  such  as  Maple  or  Mathematica.  The  Tj  matrix 
for  square  elements  is  relatively  sparse  on  all  four  edges. 

Given  any  function  that  is  expanded  in  terms  of  the  element  basis,  the  expansion  in  terms 
of  the  edge  basis  B  is  derived  as  follows: 

N 

Vi  =  'ViB=  VrTjB  =  [VjTJJ'B  s  VJjB  (17) 

j-0 


hence, 

Vi,,  =  T^V,  (18) 

After  Vi,,  has  been  computed  on  every  face  of  every  element,  the  flux  through  an  edge 
can  be  computed  without  regard  to  the  type  of  elements  or  the  orientation  of  the  coordinate 
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system  of  the  elements  that  border  the  edge.  At  an  arbitrary  edge,  illustrated  in  Figure  2,  we 
arbitrarily  designate  one  element  to  be  on  the  left  and  the  other  to  be  on  the  right.  The  two  edge 


Figure  2.  Relation  of  edge  coordinates  of  adjacent  elements. 

coordinate  systems  of  the  two  elements  always  point  in  opposite  directions.  A  function  of  the 
edge  coordinates  of  the  right  element  can  be  represented  in  terms  of  the  edge  coordinate  system 
of  the  left  element  simply  by  negating  the  odd  members  of  Vjj  .  In  the  three-dimensional  case, 
the  relation  is  not  as  trivial  because  a  rotation  may  also  be  involved. 

In  the  present  work,  the  Riemann  flux  is  approximated  by  a  simple  Lax-Friedrichs  flux  of 
the  form 

■  ds  =  +  Fr)]  -n  -  a{Vr  -  (19) 

where  ds  =  nds,  the  subscripts  /  and  r  denote  the  left  and  right  sides  of  the  edge,  n  points 
from  left  to  right,  and  a  is  some  smooth  positive  function  that  is  greater  in  magnitude  that 
the  eigenvalues  of  the  Jacobian  of  \  -I-  Ft^  ■  n  .  By  applying  equation  (18),  the 

Lax-Friedrichs  flux  is  easily  expressed  in  terms  of  the  edge  coordinate  to  give 

Ji-'^Fi^  ■  ds  =  f]b  ds  (20) 

where 

Fi  =  [/o, /,!  7,2/3 ..  .] 

=  {[jJ-i(Ty,F,  +  ITjvFr)]  -n  -  T,,V,)}/2 

=  {[■/J"^(F/  +  I^r)]  •n-a(iv,  -  Vi)}/2  (21) 
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and  I  =  diag{l^  -1,  1,  -1,  . . .)  accounts  for  the  difference  in  the  left  and  right  edge 

coordinates.  The  coefficients  of  the  approximate  Riemann  flux  on  the  right  face  are  simply 

_  _ 

Fr  =  As  a  side  note,  depending  on  the  form  of  the  flux,  less  computational  effort 

may  be  required  to  compute  F  directly  from  V  than  to  translate  the  flux  from  element  to  edge 
coordinates,  particularly  for  linear  fluxes. 

Finally,  the  boundary  integral  is  evaluated  in  the  edge  frame  of  reference  by  expressing  b}.  in 
terms  of  the  edge  coordinates  and  collecting  in  terms  of  the  components  of  F  =  [/o,  /i,  /2,  •  •  •]  • 
To  illustrate,  let  [tk]j  denote  the  fcth  row  of  Tj.  The  integral  on  a  boundary  segment  becomes 

I  ds  =  j  (lit],^B)(F’-B)<is 

dQ  dCl 

=  j  [(4,0  +  4,1^  +  4,2^^  +  4,3^  +  •••)■  (/o  +  hi  +  +  •••)] 

dQ 

—  /  [(4,0  +  4,1^  + 4,2^  +  ---)/o  +  (4,0^  +  4,1'P  +  4-,2'P  + 

dQ 

+  (4,0  ^  +  4,1^  +  4,2^^  +  -  •  •)f2  .  .  .]  (^5 

^  [ek]F  (22) 


where  [e^]  is  a  constant  row  matrix  that  is  easily  evaluated  exactly.  Let  E  denote  the  matrix 
that  is  generated  by  applying  the  above  process  to  each  member  of  the  basis  set.  The  final  form 
of  the  semidiscrete  equation  is 


V.  -  J 


-1 


M-^G 


F  (M-^Ekfk) 


k=:l 


=  0 


(23) 


where  ng  is  the  number  of  edges  and  the  matrices  M~^G  and  are  constant  matrices  that 

apply  to  all  elements  of  a  given  type.  Furthermore,  these  matrices  can  be  efficiently  precomputed 
by  the  procedure  just  described. 
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Computational  Effort 


The  effort  required  to  evaluate  the  complete  spatial  operator  is  contained  in  three  basic 
operations:  the  evaluation  of  equation  (23)  for  each  element,  the  evaluation  of  equation  (18) 
for  each  edge  of  each  element,  and  the  computation  of  F  and  F  from  V  and  V,  respectively. 
The  operation  count  of  the  first  two  operations  is  directly  related  to  the  size  of  the  matrices 
M“^G,  and  Tj.  The  row  dimension  of  all  three  matrices  is  N{p,d)  +  1.  If  the  flux 

is  linear,  then  the  column  space  of  M“^G  is  also  iV(p,  d)  +  1.  However,  in  the  nonlinear  case 
the  flux  must  be  expanded  to  at  least  degree  p  +  1 ;  thus,  the  column  dimension  must  be  at  least 
N{p  +  1,  <f)  +  1.  The  column  dimension  of  both  and  Tj  is  N{p,  d  — \  1 . 


In  a  conventional  DG  implementation  (i.e.,  one  that  uses  quadrature  points),  the  operation 
count  is  on  the  order  of 

[iV(p,d)  +  1][(1  +d)iVg^  +  2neNgb]  +  0{dNgv  +  neNgb)  (27) 

where  Ngj,  and  Ngi  denote  the  number  of  quadrature  points  required  for  the  volume  integral  and 
boundary  segment  integrals,  respectively,  and  the  last  term  denotes  the  cost  of  computing  the 
flux  at  each  quadrature  point.  Most  references  do  not  give  the  specific  quadrature  formulas  used; 
however.  Halt  et  al.^  referred  to  the  work  of  Dunavant,^^  who  derived  nearly  optimal  formulas 
in  which  Ngv  >  N{p,d)  +  I  in  order  to  evaluate  the  integral  exactly  to  degree  2p.  Thus,  the 
operation  count  for  the  conventional  DG  implementation  is  greater  than  the  values  given  by 
either  (26)  or  (27)  regardless  of  whether  we  take  advantage  of  the  sparseness  of  the  matrices. 

When  the  DG  method  is  compared  with  fundamentally  different  methods  such  as  finite- 
difference  or  finite-volume  methods,  the  comparison  must  be  done  in  an  equitable  manner.  To 
do  so,  we  hypothesize  that  any  two  methods  with  the  same  order  of  accuracy  and  the  same 
physical  stencil  size  will  yield  similar  results  (for  benign  cases  that  do  not  violate  the  basic 
assumptions  of  the  method).  In  practice,  we  compare  methods  that  are  of  the  same  order  of 
accuracy  and  have  the  same  total  number  of  dependent  variables.  In  this  frame  of  reference,  the 
evaluation  of  the  spatial  operator  is  an  operation  of  order  N{p,d)  -f  1  or  N{p  -1- 1,  d)  -f  1  per 
dependent  variable  for  a  linear  or  nonlinear  problem,  respectively. 

Time  Integration  and  Stability 

The  solution  is  advanced  in  time  with  a  three-stage  TVD  Runge-Kutta  method: 

_  yn-l 

+  (1  -  A:  =  1,2,3  (28) 

y”  = 

where  0k  =  0^  3/4,  and  1/3  for  A;  =  1,  2,  and  3,  respectively. 
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Fourier  stability  analysis  has  been  applied  to  a  generalized  K-stage  form  of  equation  (28) 


(i.e.,  k  =  1, 2,  for  the  one-dimensional  linear  case  of 


Ut  -t"  uUx  —  0  (29) 

to  determine  the  stability  limit  \j{  =  a  AtjAx.  Results  are  given  in  Table  2  for  JiT  =  1,  2,  and  3 
and  for  0  <  p  <  11.  The  rapid  drop  in  the  stability  limit  as  the  order  of  the  method  is  increased 
would  normally  be  alarming  in  comparison  with  stability  constraints  of  explicit  finite-difference 
methods.  However,  if  we  again  require  that  comparisons  be  made  among  methods  that  have  the 
same  total  number  of  variables,  then  the  size  of  the  element  in  the  DG  method  would  be  larger 
(by  a  factor  of  p  -|-  1  in  one  dimension)  than  the  mesh  size  of  a  comparable  finite-difference 
calculation.  Thus,  most  of  the  drop  in  the  stability  limit  can  be  attributed  to  the  definition  of 
Ax.  The  last  column  of  Table  2  gives  Xk{p  +  1)>  which  gives  the  DG  stability  limit  in  a  form 
that  facilitates  comparison  with  the  stability  limit  of  a  finite-difference  method. 

Results 

One-Dimensional  Test 

The  one-dimensional  version  of  this  method  has  been  tested  on  the  linear  problem 


Uf  "h  ciUx  —  0 


(30) 


and  the  nonlinear  problem 

Ut  +  liu^),  =  0  (31) 

on  the  domain  0  <  x  <  1  with  periodic  boundary  conditions. 

A  linear  problem  is  solved  first  with  smooth  initial  conditions  (7(0,  x)  =  ^  +  sin(27rx)  to 
demonstrate  the  general  accuracy  properties  of  the  method.  The  numerical  solution  is  initialized 
by  expanding  the  initial  condition  in  a  Taylor  series  about  the  center  of  each  element.  All 
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components  of  the  numerical  solution  are  compared  with  the  Taylor  series  of  the  exact  solution 
after  it  has  advected  for  several  periods.  The  L„  error  of  the  yth  component  of  the  solution  is 


defined  as 


where  uij  denotes  the  Taylor  coefficient  of  the  exact  solution  in  cell  i  and  I  denotes  the  number 
of  elements.  A  mesh-refinement  study  has  been  performed  for  p  =  1  through  5.  The  time 
step  was  chosen  to  be  sufficiently  small  such  that  the  error  would  be  dominated  by  the  spatial 
operator;  however,  for  p  >  2  the  time  step  varied  as  At  a  so  that  the  temporal 

accuracy  would  be  of  the  same  order  as  the  spatial  accuracy.  Figure  3  shows  the  Li  error  for  each 
component  of  the  solution  for  p  =  1,  2,  and  4.  The  convergence  rate  of  the  solution  between 
any  two  grids,  a  and  6,  is  defined  by 

Table  3  gives  the  convergence  rate  between  the  two  finest  grids  in  the  refinement  sequence. 
Although  most  cases  converge  at  the  design  rate  of  p  -t-  1,  the  t>o  term  of  the  p  =  1  case 
converges  at  a  rate  of  approximately  3,  which  is  1  higher  than  expected.  This  faster  convergence 
is  fortuitous  and  occurs  only  because  the  basis  functions  are  incidentally  orthogonal.  Because 
is  only  second  order  and  V2  is  undefined,  a  solution  of  degree  greater  than  1  cannot  be  recovered 
at  any  point  other  than  the  element  center  without  departing  from  the  Galerkin  framework. 
Furthermore,  although  uo  converges  faster  than  the  design  order,  its  error  is  still  considerably 
larger  than  that  of  the  formally  third-order  case  (p  =  2). 

In  the  second  test  case,  the  DG  method  is  applied  to  the  linear  problem  with  a  discontinuous 
initial  solution: 


i7(0,a;) 


4  <  ^  <  4 

®  <  i,  X  >  j 
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Figures  4  through  8  show  several  results  for  this  case  on  a  grid  with  40  elements.  Figures  4, 
5,  and  6  show  the  solution  after  one  period  for  p  =  1,  2,  and  6,  respectively.  Each  method  has 
small  overshoots  that  are  confined  to  the  neighborhood  of  the  discontinuity.  Similar  results  were 
observed  for  up  to  p  =  11  (twelth  order).  Figures  7  and  8  also  show  solutions  for  p  =  6  (seventh 
order),  but  the  solution  has  advected  for  5  and  50  periods,  respectively.  The  overshoots  neither 
grow  nor  spread  in  time,  which  is  in  sharp  contrast  to  the  behavior  of  more  traditional  methods. 
A  typical  finite-difference  approach,  for  instance,  would  tend  to  smear  a  contact  discontinuity 
over  a  region  that  grows  linearly  in  time. 

Next,  the  DG  method  is  applied  to  a  linear  test  case  that  was  prescribed  as  part  of 
the  ICASE/LaRC  Workshop  on  Benchmark  Problems  in  Computational  Aeroacoustics;*^  these 
results  are  compared  with  the  finite-difference  results  described  in  reference  14.  The  test  case 
consists  of  a  Gaussian  pulse  that  is  advected  across  a  uniform  domain.  The  Gaussian  pulse  has 
a  half-width  of  6  and  is  initially  centered  on  the  origin  of  a  domain  that  ranges  from  -20  to  450. 
Results  are  shown  in  Figure  9  for  p  =  1,2,  and  3;  however,  as  p  is  increased,  the  number  of 
elements  is  decreased  such  that  the  total  number  of  variables  is  approximately  470  (the  number 
of  points  specified  in  the  workshop).  In  Figure  10,  the  results  of  the  fourth-order  DG  method  at  t 
=  400  are  compared  in  detail  with  the  results  of  a  fourth-  and  fifth-order  finite-difference  method. 
The  results  obtained  with  the  fourth-order  DG  method  with  only  117  elements  is  considerably 
better  than  that  of  either  the  fourth-  or  the  fifth-order  finite-difference  method  with  470  points. 
(Note  that  smooth  curves  are  generated  for  the  results  obtained  with  the  DG  method  by  evaluating 
the  solution  at  several  points  within  each  element.) 

The  last  one-dimensional  test  case  is  a  nonlinear  problem  (eq.  (31))  in  which  a  shock 
forms  from  an  initially  smooth  solution.  This  problem  was  used  to  not  only  demonstrate  the 
robustness  of  the  method  but  also  to  investigate  the  effect  of  truncating  the  nonlinear  flux  at 
various  levels.  We  expect,  based  on  the  formulation,  that  the  nonlinear  flux  must  be  expanded 
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to  M  =  N{p  +  1,  c?)  terms  such  that  the  degree  of  S7bk  ■  F  is  the  same  as  the  degree  of  V  to 
obtain  the  design  rate  of  convergence  of  1.  This  expected  convergence  property  was  verified 
by  a  mesh-refinement  study  in  which  the  calculation  was  stopped  just  before  shock  formation. 
The  mesh  refinement  was  performed  for  several  values  of  p,  and  the  finest  grid  contained  320 
elements.  Typical  results,  shown  in  Table  4  for  p  =  2,  indicate  that  the  convergence  rate  of 
the  Loo  error  drops  to  approximately  p  when  M  =  N{p,  d)  but  is  approximately  p  -f  1  for 
all  other  cases. 

Figure  1 1  shows  solutions  for  p  =  1  and  p  =  2  (second  and  third  order)  in  which  the  shock 
has  formed  and  has  begun  to  propagate.  In  both  cases,  the  solutions  were  obtained  without  the 
use  of  limiters,  added  dissipation,  or  entropy  correction  terms.  However,  for  the  case  in  which 
p  =  2  the  nonlinear  flux  needed  to  be  fully  expanded  (M  =  iV(2p,  d)y,  otherwise,  the  solution 
diverges  shortly  after  shock  formation.  All  higher  order  cases  (p  >  2)  required  some  type  of 
limiter.  (Research  is  continuing  in  this  area.) 

Two-Dimensional  Test 


The  DG  method  is  applied  to  the  scalar  advection  equation  in  two  dimensions  to  demonstrate 

its  robust  treatment  of  unstructured  grids.  The  test  problem  is  given  by 

Ut  +  aU^  +  bUy  ^  f) 


(35) 

U{Q,x,y)  =  [sin(7ra:)  sin(7rp)] 

defined  on  the  periodic  domain  0  <  a:,t/  <  1.  The  approximate  solution  Vi  is  initialized  from 
the  Taylor  expansion  of  the  exact  initial  condition.  The  baseline  case  is  chosen  to  be  a  uniform 
Cartesian  grid  that  is  triangulated  in  a  regular  manner,  as  illustrated  in  Figure  12.  Figure  13 
shows  the  L\  error  in  the  uo  component  of  the  solution.  As  in  the  one-dimensional  case,  the  time 
step  was  small  so  that  the  spatial  error  dominated  and  for  p  >  2  the  time  step  was  proportional 
to  In  mesh-refinement  studies,  the  first  grid  in  the  sequence  is  coarsened  as  the 

order  of  the  method  is  increased,  so  that  the  total  number  of  variables  is  roughly  the  same. 
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The  abscissa  in  Figure  13  is  the  square  root  of  the  total  number  of  variables,  which  facilitates 
comparison  with  a  simple  fourth-order  finite-difference  method.  The  higher  order  convergence 
of  the  p  =  1  case  that  was  observed  in  one  dimension  is  not  observed  in  the  two-dimensional 
case.  The  accuracies  of  the  fourth-order  DG  and  finite-difference  methods  are  quite  similar. 

One  of  the  major  motivations  for  pursuing  a  DG  method  is  its  ability  to  maintain  accuracy 
for  complex  geometries.  Here,  the  baseline  grid  is  altered  in  several  ways  to  test  and  demonstrate 
this  capability.  Figure  14  illustrates  four  of  the  variations  that  were  tested.  In  the  first  case, 
grid  A  is  uniform  as  in  the  baseline  case,  but  the  triangulation  has  been  performed  randomly. 
In  the  second  case,  grid  B  is  generated  from  grid  A  by  smoothly  clustering  the  grid  toward  a 
diagonal.  Grid  C  is  generated  from  grid  B  by  randomly  perturbing  each  grid  point  by  an  amount 
that  is  less  than  20  percent  of  the  average  mesh  size.  In  the  last  case,  grid  D  is  generated  from 
the  baseline  by  imposing  a  piecewise-constant  mesh  spacing  that  places  half  of  the  points  in  a 
narrow  band  around  the  axis. 

In  most  cases,  the  measured  error  was  insensitive  to  the  grid  modification  or  the  direction 
of  propagation  (i.e.,  the  value  of  a  and  6).  Figure  15  gives  the  L\  error  for  p  =  3  (fourth  order) 
on  each  of  the  grids  shown  in  Figure  14.  Results  for  other  values  of  p  were  similar.  The  results 
in  Figure  15  are  for  a  =  1  and  6  =  0  (except  for  case  D2,  where  a  =  1  and  6=1).  The  slight 
increase  in  error  in  case  D2  is  attributed  more  to  the  increase  in  mesh  size  along  the  propagation 
path  than  to  the  discontinuous  manner  in  which  the  mesh  size  changes. 

In  the  last  test  case,  the  DG  method  is  applied  to  another  problem  prescribed  as  part  of  the 
ICASE/LaRC  Workshop  on  Benchmark  Problems  in  Computational  Aeroacoustics.''^  The  linear 
Euler  equations  are  solved  on  a  square  domain  of  dimensions  —100  <  x,y  <  100  with  initial 
conditions  that  place  a  compact  acoustic  source  at  x  =  y  =  0  and  a  convecting  disturbance  at 
a:  =  67,  p  =  0.  Here,  the  equations  have  been  recast  in  a  form  that  emphasizes  the  decoupling 
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of  the  convection  terms  from  the  acoustic  terms  that  occurs  in  this  linear  system: 


dE  dF 
dt  dx  dy 


(36) 


\vhere 

p-Pl  [M,ip-P)l  lMy{p-Py 

P  p,  _  MxP  +  u  p  _  MyP  +  V 

u  ’  MxU  +  P  ’  -  MyU 

V  MxV  _  MyV  +  -P  _ 

Mx  =  0.5,  My  =  0,  and 

=  0.1  exp  (-In  (2)) 

P{0,x,y)  =  exp  (-ln(2))|^^  ^ 

u{0,x,y)  =  0.04xP(0,  X,  y) 
v{0,x,y)  =  OMy  P{0,x,y) 

The  workshop,  which  targeted  finite-difference  methods,  prescribed  a  grid  of  200  x  200.  In  the 
present  calculation,  a  uniform  nxn  Cartesian  grid  is  used  that  has  been  randomly  triangulated 
(as  shown  in  Figure  14(a))  and  in  which  the  number  of  elements  is  chosen  as  a  function  of 
the  degree  p  such  that  the  total  number  of  unknowns  equals  approximately  200^.  Results  are 
shown  for  p  =  0,  1,  2,  and  3  (first,  second,  third,  and  fourth  order)  with  n  =  141,  81,  57,  and 
44,  respectively.  Figure  16  shows  P  and  u  at  i  =  40  for  the  p  =  3  case.  The  wave  fronts 
appear  smooth  and  cylindrical  in  spite  of  the  fact  that  the  initial  disturbance  was  smaller  than 
the  element  size.  A  more  quantitative  comparison  is  shown  in  Figures  17  and  18.  The  pressure 
P  is  plotted  on  the  x  =  0,  f  =  40  line  for  p  =  0  through  3  at  the  resolution  given  above. 
Also  shown  is  a  fine-grid  solution  with  p  =  3  and  n  =  132  and  the  solution  from  a  fifth-order 
finite-difference  method  on  a  200x200  grid.  An  enlargement  of  the  right  peak  (Figure  18)  shows 
that  all  solutions  of  third-order  accuracy  or  better  yield  similar  results. 


(x  -67)^ 
25 
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Conclusions 


A  quadrature-free  form  of  the  discontinuous  Galerkin  (DG)  method  has  been  formulated  for 
the  hyperbolic  conservation  laws.  This  approach  reduces  both  the  storage  and  operation  count 
to  levels  that  are  comparable  to  high-order  finite-volume  methods.  The  method  is  well  suited  to 
both  unstructured  and  structured  grids  and  has  been  tested  on  several  one-  and  two-dimensional 
problems  to  demonstrate  its  accuracy  and  robustness.  On  smooth  meshes,  the  accuracy  of  the  DG 
method  is  comparable  to  or  better  than  traditional  high-order  finite-difference  methods.  Contact 
discontinuities  are  advected  without  the  usual  diffusion  effect,  and  nonlinear  discontinuities 
(shocks)  are  propagated  by  second-  and  third-order  methods  without  the  use  of  limiters.  On 
two-dimensional  unstructured  grids,  random  and  discontinuous  mesh  variations  had  little  effect 
on  the  error  and  no  effect  on  the  convergence  of  the  error. 
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Table  1.  N(p,d)  +  1  for  Specific  Values  of  p  and  d 


d 

p 

1 

2 

3 

4 

1 

2 

3 

4 

5 

2 

3 

6 

10 

15 

3 

4 

10 

19 

31 

4 

5 

15 

31 

53 

5 

6 

21 

46 

81 

6 

7 

28 

64 

115 

Table  2.  \j{ 

III 

ft 

l> 

> 

for  ^T-stage  Runge-Kutta  Methods  Applied 

DG  Method  of  Degree  p  (order  p  +  1): 

t  denotes  unstable  method 

P 

^1 

^2 

^3 

A3(p  +  1) 

0 

1.0 

1.00 

1.256 

1.256 

1 

0.001 

0.333 

0.409 

0.818 

2 

$ 

0.06 

0.209 

0.627 

3 

$ 

0.02 

0.13 

0.52 

4 

4 

0.01 

0.089 

0.445 

5 

$ 

0.006 

0.066 

0.396 

6 

0.004 

0.051 

0.306 

7 

$ 

0.003 

0.04 

0.32 

8 

4 

0.002 

0.033 

0.297 

9 

X 

JL 

0.002 

0.027 

0.27 

10 

t 

0.001 

0.023 

0.253 

11 

t 

0.001 

0.02 

0.24 

Unstable  method. 

Table  3. 

Convergence  Rates  of  L^  (, 

Between  Two  Finest  Grids 

j 

P 

0 

1 

2  3 

4 

1 

2.990 

2.133 

2 

3.064 

3.049 

3.022 

4 

4.99 

5.74 

5.00  5.00 

5.00 
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Table  4.  Effect  of  Truncating  Nonlinear  Flux  on  Grid  Convergence:  p  =  2  case 


Convergence  Rate 

M 

N{p,d) 

N{p+l,d) 

N{p  +  2,d) 

2.894 

2.944 

2.942 

2.015 

2.923 

2.921 

(/) 

Figure  3.  Convergence  of  L\  error  for  p  =  1,  p  =  2,  and  p  =  4. 

- exact  solution 


Figure  4.  Solution  of  linear  problem  after  one  time  period,  produced  by  DG  method  with  p 
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exact  solution 
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Figure  5.  Solution  of  linear  problem  after  one  time  period,  produced  by  DG  method  with  p  =  2. 


exact  solution 
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Figure  6.  Solution  of  linear  problem  after  one  time  period,  produced  by  DG  method  with  p  =  6. 
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exact  solution 


Figure  7.  Solution  of  linear  problem  after  five  time  periods,  produced  by  DG  method  with  p  =  6. 

- exact  solution 


Figure  8.  Solution  of  linear  problem  after  50  time  periods,  produced  by  DG  method  with  p  =  6. 
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Figure  10.  Comparison  of  fourth-order  DG  method  with  fourth- 
and  fifth-order  finite-difference  (FD)  method. 
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Figure  12.  Triangulated  grid  and  solution  of  scalar  advection  problem. 
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Figure  13.  Convergence  of  solution  for  scalar  advection  on 


unstructured  grid  for  various  orders  of  accuracy. 


a)  b) 


Figure  14.  Variation  on  baseline  unstructured  grid,  a)  Grid  A,  Random  triangulation,  b)  Grid 
B,  Smoothly  clustered  toward  diagonal,  c)  Grid  C,  Random  perturbation  of  20 
percent  of  the  mean  cell  spacing,  d)  Grid  D,  Discontinuous  mesh  variation. 


Log(sqrt(/(V+l))) 

Figure  15.  Convergence  of  solution  for  scalar  advection  for  fourth-order 
method  on  variations  of  baseline  unstructured  grid. 
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